# ============================================================
# Thomas König & Stefan Eschenwecker
# The European Court of Justice and legal European integration
# 
# This script produces Figure 2 in the main text.
# ============================================================

# load packages
library(rstan)
library(tidybayes)
library(brms)
library(tidyverse)
library(patchwork)

# set plot options
theme_set(theme_minimal(base_size = 13))

# load data
CourtData <- read.csv("Court/Data/PreparedCourtData.csv")

# load models
EqVarModel <- readRDS("Court/Models/EqualPrecision.rds")
UnEqVarModel <- readRDS("Court/Models/UnequalPrecision.rds")


# prepare data
# (DV + IVs as ordered factors)
# (remove cases where all actors claimed not applicable)
CourtData <- CourtData %>%
  mutate(
    pos_cjeu_ln = factor(pos_cjeu_ln,
      ordered = T,
      levels = c("PS", "Ambi", "ME")
    ),
    supra_signal = factor(supra_signal,
      ordered = F,
      levels = c( "Uninfo", "Cont", "Weak PS", "Strong PS", "Weak ME", "Strong ME")
    ),
    supra_signal_int = relevel(factor(supra_signal_int), ref = "No Intensity")
  ) %>%
  filter(all_not_applic == 0)

# calculate randomized quantile residuals for location-only model
# (Hint: takes some time to run)
set.seed(1808)
EqResids <- CourtData %>% 
  add_predicted_draws(EqVarModel) %>%
      mutate(.prediction = ordered(levels(pos_cjeu_ln)[.prediction],
                                   levels = levels(pos_cjeu_ln))) %>%
      summarise(
        p_lower = mean(.prediction < pos_cjeu_ln),
        p_upper = mean(.prediction <= pos_cjeu_ln),
        p_residual = runif(1, p_lower, p_upper),
        z_residual = qnorm(p_residual),
        .groups = "drop_last"
)

# generate plot
EqResidsPlot <- ggplot(EqResids, aes(sample = z_residual)) +
  geom_qq(size = 1) +
  geom_abline() +
  xlab("Normal Deviates") +
  ylab("Sample Deviates")

# same for combined model
# (Hint: takes some time to run)
set.seed(1808)
UnEqResids <- CourtData %>% 
  add_predicted_draws(UnEqVarModel) %>%
      mutate(.prediction = ordered(levels(pos_cjeu_ln)[.prediction],
                                   levels = levels(pos_cjeu_ln))) %>%
      summarise(
        p_lower = mean(.prediction < pos_cjeu_ln),
        p_upper = mean(.prediction <= pos_cjeu_ln),
        p_residual = runif(1, p_lower, p_upper),
        z_residual = qnorm(p_residual),
        .groups = "drop_last"
)

UneqResidsPlot <- ggplot(UnEqResids, aes(sample = z_residual)) +
  geom_qq(size = 1) +
  geom_abline() +
  xlab("Normal Deviates") +
  ylab("Sample Deviates")


# combine plots
CombinedPlot <- EqResidsPlot +
  ggtitle("Location-Only Model") +
  UneqResidsPlot +
  ggtitle("Combined Model") +
  plot_layout(ncol = 1,
              guides = "collect",
              axis_titles = "collect")

# save plot
ggsave(plot = CombinedPlot, device = "pdf",
       filename = "ResidPlot.pdf",
       path = "Court/Plots/", height = 11, width = 8)
